inline
kmeans::kmeans(const mat &in_data, const uword in_Ncent)
: data(in_data), N_cent(in_Ncent), Dim(in_data.n_rows), N_data(data.n_cols)
{
}

inline
void
kmeans::run(const uword max_iter)
{
  err = 1;
  uword iter = 0;
  ini_cent();// Initializing cent_1
  //mat saveCent=cent_1.t();
  //saveCent.save("C_ini.dat", raw_ascii);
  cout << "running iter: ";
  while ( err > 0.005 && iter < max_iter)
  {
    update_cent(); //means are updated inside this method
    
    cout << iter << ".. ";
    
    iter++;
  }
  
  //Saving
  //saveCent=means.t();
  //saveCent.save("C_prefinal.dat", raw_ascii);
  
  //Means with less than 2 vector assigned
  uvec less2;
  less2 	= find(counts < 2); // Use only the first one
  if (less2.n_elem>0)
  { 
    //cout << "Doing check_cent";
    check_cent();//means are updated inside this method
  }
   
   
   //means.print("means in Kmeans: ");
  //Saving
  //saveCent=means.t();
  //saveCent.save("C_final.dat", raw_ascii);
  
  /*distances.zeros(N_cent);
  partitions.zeros(N_data,N_cent);
  counts.zeros(N_cent);
  
  for(uword ni = 0; ni < N_data; ++ni)
  {
    for(uword ci = 0; ci < N_cent; ++ci)
    {
      distances(ci) = norm( (data.col(ni) - means.col(ci)), 2 );
    }
    distances.min(best_ci);
    //partitions(counts(best_ci),best_ci)=ni; // Partitions is update in chec_cent
    counts(best_ci) += 1;
  }
  
  */
  //cout << "partitions after check_cent= " << endl << partitions << endl;
  //cout << "counts: "<<  counts.t() << endl;

}


inline
void 
kmeans::ini_cent()
{
  cout << "...Doing init_cent";
  uvec initial_indices;
  vec tmp = randu<vec>(N_cent) * (data.n_cols-1); 
  initial_indices = conv_to<uvec>::from( tmp );// converting to <uvec> from <vec>
  cent_1.zeros(Dim, N_cent);
  
  for(uword ci=0; ci < N_cent; ++ci)
  {
    cent_1.col(ci) = data.col( initial_indices(ci));
  }
}


inline
void 
kmeans::ini_cent(const uword Nc, const ucolvec indices) // When Cent has < 2 vectors
{
  cout << "...Doing init_cent_2";
  uvec initial_indices;
  cent_1.zeros(Dim, Nc);
  
  for(uword ci=0; ci < Nc; ++ci)
  {
    cent_1.col(ci) = data.col( indices(ci));
  }
}


inline
void 
kmeans::update_cent()
{
  
  cout << "...doing update cent" ;
  partitions.zeros(N_data,N_cent);
  distances.zeros(N_cent);
  cent_2.zeros(Dim, N_cent);
  counts.zeros(N_cent);
  
  for(uword ni = 0; ni < N_data; ++ni)
  {
    for(uword ci = 0; ci < N_cent; ++ci)
    {
      distances(ci) = norm( (data.col(ni) - cent_1.col(ci)), 2 );
    }
    distances.min(best_ci);
    partitions(counts(best_ci),best_ci)=ni;
    counts(best_ci) += 1;
  }
  
  for (uword ci=0; ci<N_cent; ++ci)
  {
    uword out=counts(ci);
    for (uword j=0; j<out ; ++j)
    {
      cent_2.col(ci)+=data.col(partitions(j,ci));
    }
    
    if (counts(ci)==0)
    {
      cent_2.col(ci).fill(0);   
    }
    else
    {
      cent_2.col(ci) /= double( counts(ci) );
    } 
    
  }
  
  double change = accu( abs( cent_1 - cent_2 ) );
  err=change;
  cent_1 = cent_2;
  means=cent_1;
}


inline
void 
kmeans::update_cent(ucolvec index) // When Cent has < 2 vectors
{
  
  //cout << "...Doing update_cent_2";
  const uword Nc = 2;
  
  
  partitions2.zeros(N_data, Nc);
  
  distances.zeros( Nc );
  
  cent_2.zeros( Dim, Nc );

  counts2.zeros(Nc);
  
  for(uword ni = 0; ni < index.n_elem; ++ni)
  {
    for(uword ci = 0; ci < Nc; ++ci)
    {
      distances(ci) = norm( (data.col(index(ni)) - cent_1.col(ci)), 2 ); // se supone que cent_1 ya se fue inicializado
    }
    distances.min( best_ci );
    partitions2( counts2(best_ci), best_ci ) = index(ni);
    counts2(best_ci) += 1;
  }
  
  //cout << "partitions2 in updating 2 Cent: "<< endl << partitions2 << endl;
  //cout << "counts2: " << counts2.t() << endl;
  for (uword ci=0; ci<Nc; ++ci)
  {
    uword out=counts2(ci);
    for (uword j=0; j<out ; ++j)
    {
      cent_2.col(ci)+=data.col(partitions2(j,ci));
    }
    
    if (counts2(ci) == 0)
    {
      cent_2.col(ci).fill(0);   
    }
    else
    {
      cent_2.col(ci) /= double( counts2(ci) );
    } 
    
  }
  double change = accu( abs( cent_1 - cent_2 ) );
  err=change;
  cent_1 = cent_2;
  //cout << "...end update_cent_2";

}


inline
void 
kmeans::check_cent()
{
  err = 1;
  uword iter = 0;
  uword maxi;
  uvec less2;
  ucolvec indices;
  less2 = find(counts < 2); // Use only the first one
  
  
   //cout << endl << endl  <<"partitions in check_cent= " << endl << partitions << endl;
   //cout << "counts: "<<  counts.t() << endl;
   //cout <<  "less2.n_elem: " << endl << less2.n_elem << endl;
  
  uvec counts_tmp=counts;
  
  for (int var=0 ; var<less2.n_elem; ++var)
  {
    counts_tmp.max(maxi);
    //cout <<  "maxi= " << maxi << endl;
    //cout << "counts_tmp(maxi)= " << counts_tmp(maxi) << endl;
    indices = partitions(span(0,counts_tmp(maxi)-1),maxi);
    //cout << "indices = " << indices.t() << endl;
    
    ini_cent(2, indices);// Initializing cent_1, with data of the most popular mean.
    while ( err > 0.005 && iter < 2)
    {
      
      update_cent(indices);
      iter++;
    }
    //cout << endl <<"2 values for cent_1.col: "  << endl;
    //cout << cent_1.col(0).t() << endl;
    //cout << cent_1.col(1).t() << endl;
    means.col(less2(var))	= cent_1.col(0);
    means.col(maxi)		= cent_1.col(1);
    //cout << "partitions:  " << endl << partitions << endl;
    //cout << "partitions2: " << endl << partitions2 << endl;
    partitions.col(less2(var))	= partitions2.col(0);
    partitions.col(maxi)	= partitions2.col(1);
    counts(less2(var)) 		= counts2(0);
    counts(maxi) 		= counts2(1);
    counts_tmp(maxi) = 0;//OJO
  }
  //cout << "end Chec_cent" << endl;
}



inline
mat
kmeans::get_means()
{
  return means;
}

inline
field <mat>
kmeans::get_covs()
{
  //cout << "kmeans::get_covs"  << endl;
  mat C;
  running_stat_vec<colvec::elem_type> stats(true);
  colvec sample;
  field<mat> covs(N_cent);
  for (uword gi=0; gi<N_cent; ++gi)
  {
    for (uword ni=0; ni<counts(gi); ++ni)
    {
      sample = data.col(partitions(ni,gi));
      //cout << "sample= " << sample << endl;
      stats(sample);
    }
    C = stats.cov();
    covs(gi) = C;
    stats.reset();
  }
  //C.diag().print("C:");
  return covs;
}




inline
vec
kmeans::get_weights()
{
  weights.zeros(N_cent);
  
  for (uword gi = 0; gi<N_cent; ++gi)
  {
    weights.at(gi) = double(counts(gi))/double(N_data);
  }
  
  //cout << "weights= " << weights.t() << endl;
  return weights;
}
